Physical properties of a GeS2 glass using approximate 
ab initio molecular dynamics 
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With the use of ab initio based molecular dynamics simula- 
tions we study the structural, dynamical and electronic prop- 
erties of glassy g-GeS2 at room temperature. From the radial 
distribution function we find nearest neighbor distances al- 
most identical to the experimental values and the static struc- 
ture factor is close to its experimental counterpart. From the 
Ge-S-Ge bond angle distribution we obtain the correct dis- 
tribution of corner and edge-sharing GeS4 tetrahedra. Con- 
cerning the dynamical characteristics we find in the mean 
square displacement of the atoms discontinuous variations 
corresponding either to the removal of coordination defects 
around a single particle or to structural rearrangements in- 
volving a larger number of atoms. Finally we calculate the vi- 
brational density of states, which exhibits two well separated 
bands as well as some features characteristic of the amorphous 
state, and the electronic density of states showing an optical 
gap of 3.27 eV. 

PACS numbers: 61.43.Fs,71.23.-k,71.15.Pd 



I. INTRODUCTION 

Amongst the chalcogenide glasses, glassy germanium 
disulfide (g-GeS2) has been heavily studied for many 
years [1] and is still the subject of recent experimen- 
tal investigations [2,3] because of its interesting physi- 
cal properties. Chalcogenide materials can be used as 
sensitive media for optical recording, as light guides, as 
high-resolution inorganic photoresistors or anti reflection 
coatings [4] . Moreover bulk glasses with for example Ag+ 
cations are good solid electrolytes with a high ionic con- 
ductivity at room temperature [5] and thin GeS2 films 
are promising materials for submicron lithography when 
doped with silver [6]. The high quantum efficiency of 
these glasses appears as a consequence of the relative 
high masses of the elements involved [7]. All these po- 
tential applications of glassy GeS2 have led many authors 
to study the physical properties of these chalcogenide 
glasses and many experiments have been done on this 
topic [8,9]. However in order to understand the physical 
mechanisms occurring at the atomic scale and leading 
to the results observed in experiments, numerical simu- 
lations can be an alternative tool and more specifically 
Molecular Dynamics (MD) simulations. Although clus- 
ter modeling simulations were performed on g-GeS2 [10], 



it appears that GeS2 compounds have not been the topic 
of extensive MD investigations yet, contrary to GeSe2 
[11] or SiSe2 [12]. In order to perform such investigations 
one has to decide what kind of description (classical or ab 
initio) is adequate for GeS2. Taking the mostly (but not 
purely) covalent bonding into account in g-GeS2 a first- 
principles approach seems appropriate. In this paper we 
present therefore a theoretical study of the structural, 
dynamical and electronic properties of g-GeS2 using an 
approximate ab initio description based on the Sankey- 
Niklewski scheme [13] and materialized in the so-called 
"FIREBALL96" MD code [14]. This technique has been 
successfully used in the study of several different chalco- 
genide glasses [15,16] and in order to check its vahdity 
in the case of GeS2 samples we have compared our re- 
sults with experimental results when those were available. 
Concerning the structure at 300K, the nearest-neighbor 
distances as well as the static structure factor compare 
well with the experimental data. Using the angle dis- 
tributions and the radial pair distribution functions we 
find the correct proportion of edge and corner sharing 
GeS4 tetrahedra which are the basic building blocks of 
the germanium disulfide glass. Concerning the dynamics 
of the individual particles, we find in the mean square 
displacement (MSD) signatures of individual or collec- 
tive atomic rearrangements corresponding to either the 
removal of "defects" or to the oscillation of large clusters 
which could be at the origin of the excess of modes seen 
at low frequency in the vibrational spectrum. 
The paper is organized as follows: In section II we present 
the theoretical foundation of the FIREBALL96 code as 
well as the approximations used. In section III the results 
are presented for the structural, dynamical and electronic 
properties of the GeS2 sample and section IV gives the 
major conclusions. 



II. MODEL 

The theoretical framework of our work is the widely 
used Density Functional Theory (DFT) [17] using three 
additional approximations. 

First, we use the well known Local Density Approx- 
imation (LDA) [18] combined with the pseudopoten- 
tial approximation, which replaces the core electrons 
by an effective potential acting on the valence elec- 
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trons (Hamman-Schluter-Chiang pseudopotentials are 
used [19]). The electronic eigenstates are determined by 
a tight-binding-like linear combination of Pseudo Atomic 
Orbitals (PAO's) that satisfy the atomic self-consistent 
Hohenberg-Kohn-Sham equations [20]. A minimal basis 
set of one s and three p confined pseudo-orbitals per site 
is required. 

The second approximation has been suggested by Har- 
ris [21]. It consists in using a sum of neutral- atom 
spherical charge densities as a 0*''-order approximation 
to the self-consistent density, keeping only the first-order 
changes from this density in the energy functional. This 
approximation avoids the necessity of iterating to self- 
consistency, so eigenvalues only need to be determined 
once instead of «10 times at each step. This approxi- 
mation also avoids four-center Coulomb integrals in our 
calculations, which is a great simplification. The Har- 
ris functional has been used in many studies and has 
always given surprisingly good agreement with fully self- 
consistent calculations, except for highly ionic systems 
[15,22]. 

A third approximation is made to reduce the range of 

the tight-binding- like Hamiltonian matrix elements. To 
that purpose, the PAO's are slightly excited by imposing 
the boundary condition that they vanish outside a pre- 
determined radius. This cut-off radius is chosen equal 
to 5ao, which represents a distance of 2.645 A. Atoms 
do not overlap each other beyond twice this distance, so 
the number of neighbors of each atom is considerably re- 
duced. 

All these approximations permit to gain a considerable 
amount of CPU time compared to ab initio methods 
like the Car-Parrinello scheme [23] and therefore one 
can perform longer simulation runs or study larger sys- 
tems. Moreover this method has proved to be a very 
efficient tool for a wide variety of problems, and has 
been used with success in many different investigations 
[15,16,24,25]. 

Concerning the details of the present simulation, all of 
our calculations were performed in the microcanonical 
ensemble (N,V,E=constant), with a time step of 2.5 fs, 
and using only the F point to sample the Brillouin zone. 
The initial configuration of our system was a crystalline 
Q;-GeS2 sample containing 96 particles (32 Ge and 64 
S) confined in a cubic cell of 13.82 A to which periodic 
boundary conditions have been applied. This represents 
a density of 2.75 g.cm^^ , which is the usual experimen- 
tal density [9]. This crystalline configuration was then 
melted at 2000 K over approximately 2 ps and then equi- 
librated at 1000 K for an additional 1.5 ps. We then 
quenched the system (by velocity rescaling) through the 
glass transition (Tg=710 K) to a target temperature of 
T=300K over 4ps (for more details on similar systems, 
see [16]). Starting from this configuration, we performed 
a very long thermal MD simulation at 300 K over 375 ps 
i.e. 150000 steps. During this time, we saved the config- 
urations every 20 steps and consequently all the results 
presented below have been averaged over these 7500 con- 



figurations. 



III. RESULTS 

A. Structural properties 

The basic building blocks of glassy g-GeS2 are GeS4 
tetrahedra, connected together within a random network. 
The structural unit disorder is refiected in the absence of 
long range order and in the wide distribution of bond 
lengths and bond angles. Structural information may be 
extracted from the radial pair correlation function g{r). 
For a given a — P pair it is defined by: 



9{r)a-f3 = 



V 



Anr'^Nadr 



(1) 



Results are shown in Fig. 1 for the three different 
pairs. The smallest distance appears for the Ge-S pairs 
(Fig. 1(b)) at 2.22 A, and is in perfect agreement with 
the distance determined experimentally (2.21 A [8]). The 
distance between two Ge atoms represents the intertctra- 
hedral distance, and depends on the nature of the connec- 
tion between the tetrahedra. The first peak at 2.91 A in 
Fig. 1(a) is due to edge-sharing tetrahedra, while the sec- 
ond one, at 3.41 A, is due to corner-sharing links as shown 
in Fig. 2. The experimental distances are respectively es- 
timated at 2.91 A and 3.42 A [8], which is extremely close 
to our results. Finally the S-S pairs are responsible of the 
wide peak centered at 3.66 A (Fig. 1(c)), which is also 
extremely close to the experimental first S-S distance of 
3.64 A [8]. 

A complementary way to analyze the structure is to com- 
pute the static structure factor S{q) (obtained by Fourier 
transformation of gir)) which can be directly compared 
to its experimental counterpart. In Fig. 3 we present the 
calculated S{q) together with the one obtained by neu- 
tron diffraction experiments [26]. The good agreement 
between the two curves shows the quality of the model 
concerning the structural description of GeS2 glasses. 
The First Sharp Diffraction Peak (FSDP), which is a 
signature of the intermediate range order in amorphous 
states, appears at wl A~^ and is slightly underestimated 
compared to the experimental one. This is probably a 
consequence of the small size of our system: 1 A~^ rep- 
resents in real space a distance of 6.3 A and a sphere with 
such a radius provides a volume which is close to the to- 
tal volume of our cell. Therefore the lack of statistics for 
these large distances can explain the underestimation of 
the FSDP in our simulation. 

In order to analyze completely the medium-range struc- 
ture we have also calculated the bond angle distributions 
and in particular the intratetrahedral SGeS and interte- 
trahedral GeSGe bond angles which are represented in 
Fig. 4. The intratetrahedral angle SGeS is centered at 
110°, which is close to the perfect tetrahedral angle of 
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109.47°. Its large distribution is a signature of the struc- 
tural disorder of our glassy sample. The intertetrahedral 
bond angle GeSGe is the angle between tetrahedra and 
inchidcs two major contributions. The first one, cen- 
tered at 80°, is caused by edge-sharing tetrahedra. The 
second, at RilOO°, is due to corner-sharing tetrahedra. 
The integration of these two peaks permits to estimate 
the fraction of edge-sharing and corner-sharing connec- 
tions which are respectively 18.6% and 81.4%. These 
results have also been confirmed by a direct counting of 
each type of connection in our sample. Experimental Ra- 
man scattering measurements in amorphous g-GcS2 have 
given 16.6% of edge-sharing links and 83.4% of corner- 
sharing links [3] which is relatively close to our results. 
In view of all these data, we can safely say that the model 
describes correctly the structure of amorphous GeS2. It 
remains to be seen if this is also true for the dynamical 
properties which is the topic of the next section. 

B. Dynamical properties 

The dynamical properties of glassy g-GeS2 have been 
studied through the Mean Square Displacement (MSD) 
and the vibrational density of states (VDOS). The MSD 
is defined as — Unit) - fl{0)\^) where fl{t) is 

the position of particle i at time t. We can deduce from 
the slope of the MSD the diffusion constant D, where 
D = i limj^cc '' I*"* • In our calculation, D was found 
equal to zero; this means that the thermal energy at am- 
bient temperature is not high enough to reach the diffu- 
sive regime during the time scale of our simulation (375 
ps). Nevertheless during this time some specific struc- 
tural rearrangements can occur which manifest them- 
selves by a brutal increase of the total MSD or of the 
MSD of individual atoms. In this later case the "jumps" 
in the MSD are due to the removal of a coordination 
"defect" in the glassy system. An example of such a re- 
arrangement is shown in Fig. 5: Fig. 5(a) represents the 
individual MSD of Gc*, a particular germanium atom, 
with a brutal increase from « 0.5 (before the jump) to 
« 2.5 (after the jump) around 100 ps. The jump can 
clearly be seen in Fig. 5(b) which shows the projection of 
the displacement of Ge* on the x-z plane while the rea- 
son of the jump becomes apparent in Fig. 5(c) and 5(d) 
which illustrates the group of particles surrounding Ge* 
just before and after the jump. Indeed we see that the 
initially 3-coordinated Ge* atom gets linked with a termi- 
nal sulfur atom creating thus 2 edge-sharing tetrahedra 
which is a configuration energetically more favorable. 
The second kind of rearrangement illustrated in Fig. 6 
involves a larger number of particles and manifests it- 
self by a "pulse" in the total MSD whose amplitude is 
more important for the sulfur atoms than for the ger- 
manium atoms as shown in Fig. 6(a). In that case a 
group of particles (w 20) in a certain configuration at ti 
(Fig. 6(b)) switches to a new state at t2 (Fig. 6(c)) which 



can be called metastable since its life time is relatively 
short (« 10 ps) before the system comes back again to 
its original structure. Note that in this case no link has 
been broken or created. We observed such "oscillations" 
twice in our simulation with a time interval of 300 ps. 
Clearly our simulation time is too short to see if these 
oscillations repeat themselves at a well determined (low) 
frequency and to make a connection with the so-called 
"soft"-modes [28] well known in amorphous systems. 
To complete the study of the dynamical properties we 
have computed giv), the Vibrational Density Of States 
(VDOS), via a Fourier transformation of the velocity au- 
tocorrelation function: 

1 f°° 

g{v) = ——'Einii / exp{ivt){vi{t).Vi{Q))dt (2) 

The Fourier transformation has been calculated using 
the Wiener-Kinchin theorem [27] over the last 4096 steps 
of the simulation. The total spectrum as well as the par- 
tial contributions due to Ge and S are shown in Fig. 7. 
Despite serious efforts we could not find the experimen- 
tal counterpart of the total spectrum since apparently 
no neutron diff'raction studies have been performed on 
g-GcS2. But comparing our results with those obtained 
for the analogous GeSe2 glasses [16] the spectrum ex- 
hibits the same features. Mainly two bands can be dis- 
tinguished: a low-energy acoustic band involving mainly 
extended interblock vibrations and a high-energy optic 
band consisting of more localized intrablock vibrations. 
The two main bands are clearly separated and have ap- 
proximately the same width (7 Thz). 

In addition to the usual acoustic and optical bands, a 
small band can be seen close to 8 Thz corresponding to 
the so-called Ai mode [15]. The Ai mode is well known 
to be a tetrahcdral breathing mode (in which a central Ge 
atom is stationary and its four S neighbors move radially 
relative to the fixed Ge). This feature is strongly revealed 
in Raman measurements [29], because the mode is espe- 
cially Raman active. In Raman measiirements, there is a 
clear indication of a "two peak" structure to the Ai band. 
In particular, one usually sees a high frequency peak or 
shoulder which is interpreted as arising from edge-sharing 
tetrahedra (see Fig. 2), and the main band from tetra- 
hedra in corner-sharing conformations [30]. The Ai and 
Ale niodes have also been resolved in inelastic neutron 
scattering studies of g-GeSe2 [31]. In our work the tetra- 
hedral breathing band does show a clearly resolved split- 
ting. It is possible that a direct analysis of the eigen- 
vectors of the dynamical matrix [31] would provide more 
information linking the observed spectral feature to mi- 
croscopic vibrational excitations. At the low frequency 
end of the spectrum, a shoulder is present between 1 and 
2 THz which is coherent with the existence of a "Boson" 
peak foimd experimentally [32]. The "Boson" peak refers 
to an excess in the VDOS with respect to the Debye dis- 
tribution and is located generally around 1.5THz. 
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C. Electronic properties 

In Fig. 8 we present the electronic density of states 
(EDOS), obtained by binning the density functional elec- 
tron energy eigenvalues from the starting, fully relaxed 
model. The F point optical gap is 3.27 eV which com- 
pares very well with the experimental value of 3.2 eV ob- 
tained by resonant Raman-scattering spectroscopy [32]. 
This good agreement is due to the opposite effects of the 
use of a minimal basis set which is well known to exag- 
gerate the gap and of Kohn-Sham eigenvalues which tend 
to underestimate the gap. It should be noted that this 
gap is greater than that (2.3 eV) of g-GeSe2. Another 
point from Fig. 8 is the lack of any localized states in 
the optical gap (the Fermi level is at i? = in our cal- 
culation). This lack of gap states is realistic, since the 
density of gap states is very small in g-GeS2. 

IV. CONCLUSION 

We find that the results obtained for our 96-atom GeS2 
model are in excellent agreement with all the correspond- 
ing experimental results that are available. This real- 
ism is surprising since the size of our system is rela- 
tively small, and accurate MD simulations usually require 
larger systems. 

The structural properties of g-GeS2, which have been ex- 
tensively studied, are all extremely realistic in our simu- 
lation. The pair-correlation functions lead to interatomic 
distances that are within 10~^A compared to the experi- 
mental values and the static structure factor is very sim- 
ilar to the one obtained from neutron diffraction studies. 
The small underestimation of the FSDP encourages us to 
use larger models and we are currently preparing samples 
containing 258 atoms. The fraction of edge and corner- 
sharing tetrahedra, which can be deduced from the angu- 
lar distribution, is also close to experiment. It should be 
mentioned that we don't find homopolar (Ge-Ge or S-S) 
bonds in the present investigation but their existence can 
not be excluded a priori in the 258-atom model. Proba- 
bly a more in-depth study of the large system will permit 
to solve the apparent disagreement between two recent 
experimental studies on this point [2,3]. 
Concerning the dynamical properties of our sample we 
find discontinuous atomic displacements at ambient tem- 
perature, leading to jumps in the MSD. These jumps 
can either be due to the removal of coordination defects 
around a single atom or to oscillations of larger groups of 
atoms (« 20) between a stable and metastable configu- 
ration which could be at the origin of "soft" -modes that 
are often seen in amorphous systems. The vibrational 
density of states of glassy GeS2 could not be compared 
directly to the experimental spectrum since to our knowl- 
edge it is not available in the literature yet. We find ba- 
sically two bands separated by a "gap" in which exists a 
small structure due to the tetrahedral breathing modes. 



At low frequency we find at around 1 Thz a shoulder cor- 
responding to the famous Boson peak present in many 
amorphous systems. Concerning the electronic proper- 
ties we find an optical gap of 3.27 eV and no localized 
states in the gap which is in agreement with experimen- 
tal data. 

These first MD results obtained for g-Ge2 with the use of 
the Fireball96 code show that one can have confidence in 
this "pscudo" ah initio scheme giving an excellent de- 
scription of the physical characteristics of germanium 
disulfide for a relatively low computer load. 
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FIG. 1. Radial pair distribution functions: (a) Ge-Ge, 
(b) Ge-S, (c) S-S 
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Fig 5 




FIG. 4. Bond angle distributions. ^^G. 5. (a) Mean Square Displacement of Ge* before and 

after the jump; (b) Projection on the (x,z) plane of the tra- 
jectory of Ge*; Atomic configuration around Ge* before (c) 
and after (d) the jump. 
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Fig 8 
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FIG. 8. Electronic density of states. Dashed vertical line is 
the Fermi energy. The optical gap is found to be 3.27 eV. 
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